Myeloid deconvolution assessment
Marc Elosua-Bayes
abril 25, 2022
Introduction
In this script we are going to map Myeloid subtypes onto the Visium slides.
Libraries
library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)
library(SPOTlight)
library(SPATA2)Setting parameters
Loading necessary paths and parameters
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
"{myeloid}/{plt_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
"{myeloid}/{robj_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
set.seed(123)Extract sample id and get Donor ID
# sample_id <- params$sample_id
sample_id <- "esvq52_nluss5"
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>% dplyr::pull(donor_id)Load data
We have 8 different datasets that we are going to analyze separately. The spatial data comes from the script 03-clustering/03-clustering_integration.Rmd while the sc data can be found in Ramon’s scRNAseq analysis: /scratch/devel/rmassoni/tonsil_atlas_private/2-DOWNSTREAM_PROCESSING/results/R_objects/processed_seurat_objects/processed_seurat_objects/tonsil_integrated_with_harmony_scrublet_annotated.rds.
sp_obj <- "misc/{robj_dir}/20220215_tonsil_atlas_spatial_seurat_obj.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# Load SPOTlight data
spotlight_ls <- "{myeloid}/{robj_dir}/spotlight_ls_myeloid.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# Load names df data
nm_df <- "{myeloid}/{robj_dir}/myeloid_nm_df.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)Add colors to cell types
# (nm_df <- data.frame(col_vec))
# nm_df$plt_nm <- rownames(nm_df)
# nm_df$mod_nm <- stringr::str_replace_all(
# string = rownames(nm_df),
# pattern = "[[:punct:]]|[[:space:]]",
# replacement = ".")Analysis
Preprocess data
decon_mtrx <- spotlight_ls[[2]]
decon_mtrx <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
# Set as 0 cell types predicted to be under 3 % of the spot
decon_mtrx[decon_mtrx < 0.03] <- 0Change column names
new_cn <- data.frame(mod_nm = colnames(decon_mtrx)) %>%
dplyr::left_join(nm_df, by = "mod_nm") %>%
dplyr::mutate(plt_nm = dplyr::if_else(is.na(plt_nm), mod_nm, plt_nm)) %>%
dplyr::distinct() %>%
dplyr::pull(plt_nm)
colnames(decon_mtrx) <- new_cnWe are going to add the deconvolution to the Seurat object.
sp_obj@meta.data <- cbind(sp_obj@meta.data, decon_mtrx)Subset sample of interest
sp_sub <- sp_obj[, sp_obj@meta.data$gem_id == sample_id]
sp_sub@images <- sp_sub@images[Seurat::Images(sp_sub) == sample_id]Look at SPOTlight results
Check Topic profiles
nmf_mod_ls <- spotlight_ls[[1]]
nmf_mod <- nmf_mod_ls[[1]]
h <- NMF::coef(nmf_mod)
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
h = h,
train_cell_clust = nmf_mod_ls[[2]])
topic_profile_plts[[2]] +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90),
axis.text = ggplot2::element_text(size = 12))
Look at all cells profiles
topic_profile_plts[[1]] +
ggplot2::theme(
axis.text.y = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank(),
axis.title = ggplot2::element_blank())
Look at cells topic profile
basis_spotlight <- data.frame(NMF::basis(spotlight_ls[[1]][[1]]))
train_labs <- spotlight_ls[[1]][[2]]
colnames(basis_spotlight) <- unique(stringr::str_wrap(train_labs, width = 30))
basis_spotlight[basis_spotlight < 0.0000001] <- 0
DT::datatable(basis_spotlight, filter = "top")Cell type location
Look at the location of each cell type in each slice separately
# Iterate over cell types
ct <- colnames(decon_mtrx)
# Iterate over images
lapply(names(sp_obj@images), function(nm) {
print(nm)
nm_donor <- id_sp_df %>% dplyr::filter(gem_id == nm) %>% dplyr::pull(donor_id)
# Iterate over cell types
ct_plt_ls <- lapply(ct, function(i) {
tmp_plt <- Seurat::SpatialFeaturePlot(
object = sp_obj,
features = i,
alpha = c(0, 1),
images = nm) +
ggplot2::scale_fill_gradientn(
colors = heat.colors(10, rev = TRUE)) +
ggplot2::scale_alpha(range = c(0, 1)) +
ggplot2::labs(title = stringr::str_wrap(string = i,
width = 25),
fill = "") +
ggplot2::theme(
plot.title = ggplot2::element_text(
hjust = 0.5,
size = 20,
face = "bold"))
if (sum(sp_sub@meta.data[, i]) == 0) {
tmp_plt <- suppressMessages(tmp_plt + ggplot2::scale_alpha(range = c(0,0)))
}
return(tmp_plt)
})
(plt_arr <- cowplot::plot_grid(
plotlist = ct_plt_ls,
axis = "trbl",
align = "hv",
nrow = 6,
ncol = 5))
"{myeloid}/{plt_dir}/cell_type_location_myeloid_{nm_donor}_new.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = plt_arr,
base_height = 30,
base_width = 25)
})## [1] "tarwe1_xott6q"## [1] "c28w2r_7jne4i"## [1] "esvq52_nluss5"## [1] "p7hv1g_tjgmyj"## [1] "gcyl7c_cec61b"## [1] "zrt7gl_lhyyar"## [1] "qvwc8t_2vsr67"## [1] "exvyh1_66caqq"## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-2-T_new.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-8-T_new.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-10-T_new.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-12-T_new.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-13-T_new.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-14-T_new.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-9-T_new.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-11-T_new.pdf"Boxplots
Prepare data for boxplots
metadata_long <- sp_sub@meta.data %>%
# tidyr::pivot_longer(cols = c("annotation"),
# names_to = "stratification_id",
# values_to = "stratification_val") %>%
tidyr::pivot_longer(cols = all_of(ct), names_to = "ct_key", values_to = "ct_val") %>%
# dplyr::left_join(col_df, by = c("ct_key" = "ct_name")) %>%
dplyr::mutate(ct_val = dplyr::if_else(ct_val > 0.001, ct_val, 0))Box plot of cell type proportion between stratified regions showing the unadjusted ANOVA Pvalue
keep_ct <- metadata_long %>%
dplyr::group_by(ct_key) %>%
dplyr::summarise(prop_sum = sum(ct_val)) %>%
dplyr::filter(prop_sum > 0) %>%
dplyr::pull(ct_key)
(bplt <- metadata_long %>%
dplyr::filter(ct_key %in% keep_ct) %>%
dplyr::mutate(
ct_key = stringr::str_wrap(string = ct_key,
width = 30)) %>%
ggpubr::ggboxplot(data = .,
x = "annotation_20220215",
y = "ct_val",
facet.by = "ct_key",
color = "annotation_20220215",
fill = "annotation_20220215",
add = "jitter",
scales = "free",
repel = TRUE,
outlier.shape = NA,
alpha = 0.6,
palette = "Set1",
ncol = 6) +
ggplot2::theme(
strip.text = ggplot2::element_text(size = 18, face = "bold"),
axis.text.y = ggplot2::element_text(size = 16),
axis.title.y = ggplot2::element_text(size = 20),
# axis.text.x = element_text(size = 12, angle = 90, vjust = 0.5, hjust = 0.5),
axis.text.x = ggplot2::element_blank(),
legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_blank(),
strip.background = ggplot2::element_blank()) +
ggplot2::labs(
y = "Proportion",
color = "Regions",
fill = "Regions"))
# bplt <- bplt +
# ggpubr::stat_compare_means(method = "anova", size = 6) +
# ggplot2::scale_y_continuous(
# expand = expansion(mult = c(0, 0.1)),
# labels = function(x) sprintf("%.2f", x))
"{myeloid}/{plt_dir}/strat_bplot_{donor_id}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = bplt,
base_height = 20,
base_width = 25)Cell type correlation matrix
Slide-specific
lapply(Images(sp_obj), function(i) {
se_sub <- subset(sp_obj, subset = gem_id == i)
# se_sub
se_sub@images <- se_sub@images[Seurat::Images(se_sub) == i]
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>%
dplyr::pull(donor_id)
(cor_mtrx_ct <- SCrafty::correlation_heatmap(
se = sp_obj,
feats = ct,
assay = "Spatial",
slot = "data") +
ggplot2::labs(
title = glue::glue("{donor_id}: Cell-type correlation matrix")))
"{myeloid}/{plt_dir}/cor-mtrx_cell-type_{donor_id}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_mtrx_ct,
base_height = 9,
base_width = 10)
})## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"Integrated
# se_sub <- subset(merged_se, subset = gem_id == "esvq52_nluss5")
# se_sub
# se_sub@images <- se_sub@images[Seurat::Images(se_sub) == "esvq52_nluss5"]
(cor_mtrx_ct <- SCrafty::correlation_heatmap(
se = sp_obj,
feats = ct,
assay = "Spatial",
slot = "data") +
ggplot2::labs(
title = "Integrated cell-type correlation matrix"))
"{myeloid}/{plt_dir}/cor-mtrx_cell-type_integrated.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_mtrx_ct,
base_height = 9,
base_width = 10)Session Info
sessionInfo()## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
##
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/singlecell/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_1.2.0 tibble_3.1.6 Biobase_2.50.0 BiocGenerics_0.36.1 SPATA2_0.1.0 SPOTlight_0.1.7 readr_2.1.2 stringr_1.4.0 dplyr_1.0.8 cowplot_1.1.1 ggpubr_0.4.0 ggplot2_3.3.5 SeuratObject_4.0.4 Seurat_4.1.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.24 tidyselect_1.1.2 htmlwidgets_1.5.4 grid_4.0.1 Rtsne_0.15 munsell_0.5.0 codetools_0.2-18 ragg_1.2.2 ica_1.0-2 DT_0.21 future_1.24.0 miniUI_0.1.1.1 withr_2.5.0 spatstat.random_2.1-0 colorspace_2.0-3 highr_0.9 knitr_1.37 stats4_4.0.1 SingleCellExperiment_1.12.0 ROCR_1.0-11 ggsignif_0.6.3 tensor_1.5 listenv_0.8.0 NMF_0.23.0 MatrixGenerics_1.2.1 labeling_0.4.2 GenomeInfoDbData_1.2.4 polyclip_1.10-0 bit64_4.0.5 farver_2.1.0 rprojroot_2.0.2 parallelly_1.30.0 vctrs_0.3.8 generics_0.1.2 xfun_0.30 R6_2.5.1 doParallel_1.0.17 GenomeInfoDb_1.26.7 bitops_1.0-7 spatstat.utils_2.3-0 DelayedArray_0.16.3
## [43] assertthat_0.2.1 promises_1.2.0.1 scales_1.1.1 vroom_1.5.7 gtable_0.3.0 globals_0.14.0 goftest_1.2-3 rlang_1.0.2 systemfonts_1.0.4 splines_4.0.1 rstatix_0.7.0 lazyeval_0.2.2 spatstat.geom_2.3-2 broom_0.7.12 yaml_2.3.5 reshape2_1.4.4 abind_1.4-5 crosstalk_1.2.0 backports_1.4.1 httpuv_1.6.5 tools_4.0.1 gridBase_0.4-7 ellipsis_0.3.2 spatstat.core_2.4-0 jquerylib_0.1.4 RColorBrewer_1.1-2 ggridges_0.5.3 Rcpp_1.0.8 plyr_1.8.6 zlibbioc_1.36.0 purrr_0.3.4 RCurl_1.98-1.6 rpart_4.1.16 deldir_1.0-6 pbapply_1.5-0 S4Vectors_0.28.1 zoo_1.8-9 SummarizedExperiment_1.20.0 ggrepel_0.9.1 cluster_2.1.2 here_1.0.1 magrittr_2.0.2
## [85] data.table_1.14.2 scattermore_0.8 lmtest_0.9-39 RANN_2.6.1 fitdistrplus_1.1-6 matrixStats_0.61.0 hms_1.1.1 patchwork_1.1.1 mime_0.12 evaluate_0.15 xtable_1.8-4 IRanges_2.24.1 gridExtra_2.3 compiler_4.0.1 KernSmooth_2.23-20 crayon_1.5.0 htmltools_0.5.2 mgcv_1.8-39 later_1.3.0 tzdb_0.2.0 DBI_1.1.2 corrplot_0.92 MASS_7.3-55 Matrix_1.4-0 car_3.0-12 cli_3.2.0 SCrafty_0.1.0 igraph_1.2.11 GenomicRanges_1.42.0 pkgconfig_2.0.3 registry_0.5-1 plotly_4.10.0 spatstat.sparse_2.1-0 foreach_1.5.2 ggcorrplot_0.1.3 bslib_0.3.1 rngtools_1.5.2 pkgmaker_0.32.2 XVector_0.30.0 digest_0.6.29 sctransform_0.3.3 RcppAnnoy_0.0.19
## [127] spatstat.data_2.1-2 rmarkdown_2.12 leiden_0.3.9 uwot_0.1.11 shiny_1.7.1 lifecycle_1.0.1 nlme_3.1-155 jsonlite_1.8.0 carData_3.0-5 viridisLite_0.4.0 fansi_1.0.2 pillar_1.7.0 lattice_0.20-45 fastmap_1.1.0 httr_1.4.2 survival_3.3-1 glue_1.6.2 png_0.1-7 iterators_1.0.14 bit_4.0.4 stringi_1.7.6 sass_0.4.0 textshaping_0.3.6 irlba_2.3.5 future.apply_1.8.1